{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Source of the materials**: Biopython cookbook (adapted)\n",
    "<font color='red'>Status: Draft</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BLAST and other sequence search tools (*experimental code*) {#chapter:searchio}\n",
    "===========================================================\n",
    "\n",
    "*WARNING*: This chapter of the Tutorial describes an *experimental*\n",
    "module in Biopython. It is being included in Biopython and documented\n",
    "here in the tutorial in a pre-final state to allow a period of feedback\n",
    "and refinement before we declare it stable. Until then the details will\n",
    "probably change, and any scripts using the current `Bio.SearchIO` would\n",
    "need to be updated. Please keep this in mind! For stable code working\n",
    "with NCBI BLAST, please continue to use Bio.Blast described in the\n",
    "preceding Chapter \\[chapter:blast\\].\n",
    "\n",
    "Biological sequence identification is an integral part of\n",
    "bioinformatics. Several tools are available for this, each with their\n",
    "own algorithms and approaches, such as BLAST (arguably the most\n",
    "popular), FASTA, HMMER, and many more. In general, these tools usually\n",
    "use your sequence to search a database of potential matches. With the\n",
    "growing number of known sequences (hence the growing number of potential\n",
    "matches), interpreting the results becomes increasingly hard as there\n",
    "could be hundreds or even thousands of potential matches. Naturally,\n",
    "manual interpretation of these searches’ results is out of the question.\n",
    "Moreover, you often need to work with several sequence search tools,\n",
    "each with its own statistics, conventions, and output format. Imagine\n",
    "how daunting it would be when you need to work with multiple sequences\n",
    "using multiple search tools.\n",
    "\n",
    "We know this too well ourselves, which is why we created the\n",
    "`Bio.SearchIO` submodule in Biopython. `Bio.SearchIO` allows you to\n",
    "extract information from your search results in a convenient way, while\n",
    "also dealing with the different standards and conventions used by\n",
    "different search tools. The name `SearchIO` is a homage to BioPerl’s\n",
    "module of the same name.\n",
    "\n",
    "In this chapter, we’ll go through the main features of `Bio.SearchIO` to\n",
    "show what it can do for you. We’ll use two popular search tools along\n",
    "the way: BLAST and BLAT. They are used merely for illustrative purposes,\n",
    "and you should be able to adapt the workflow to any other search tools\n",
    "supported by `Bio.SearchIO` in a breeze. You’re very welcome to follow\n",
    "along with the search output files we’ll be using. The BLAST output file\n",
    "can be downloaded\n",
    "[here](http://biopython.org/SRC/Tests/Tutorial/my_blast.xml), and the\n",
    "BLAT output file\n",
    "[here](http://biopython.org/SRC/Tests/Tutorial/my_blat.psl). Both output\n",
    "files were generated using this sequence:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    ">mystery_seq\n",
    "CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The BLAST result is an XML file generated using `blastn` against the\n",
    "NCBI `refseq_rna` database. For BLAT, the sequence database was the\n",
    "February 2009 `hg19` human genome draft and the output format is PSL.\n",
    "\n",
    "We’ll start from an introduction to the `Bio.SearchIO` object model. The\n",
    "model is the representation of your search results, thus it is core to\n",
    "`Bio.SearchIO` itself. After that, we’ll check out the main functions in\n",
    "`Bio.SearchIO` that you may often use.\n",
    "\n",
    "Now that we’re all set, let’s go to the first step: introducing the core\n",
    "object model.\n",
    "\n",
    "The SearchIO object model {#sec:searchio-model}\n",
    "-------------------------\n",
    "\n",
    "Despite the wildly differing output styles among many sequence search\n",
    "tools, it turns out that their underlying concept is similar:\n",
    "\n",
    "-   The output file may contain results from one or more search queries.\n",
    "\n",
    "-   In each search query, you will see one or more hits from the given\n",
    "    search database.\n",
    "\n",
    "-   In each database hit, you will see one or more regions containing\n",
    "    the actual sequence alignment between your query sequence and the\n",
    "    database sequence.\n",
    "\n",
    "-   Some programs like BLAT or Exonerate may further split these regions\n",
    "    into several alignment fragments (or blocks in BLAT and possibly\n",
    "    exons in exonerate). This is not something you always see, as\n",
    "    programs like BLAST and HMMER do not do this.\n",
    "\n",
    "Realizing this generality, we decided use it as base for creating the\n",
    "`Bio.SearchIO` object model. The object model consists of a nested\n",
    "hierarchy of Python objects, each one representing one concept outlined\n",
    "above. These objects are:\n",
    "\n",
    "-   `QueryResult`, to represent a single search query.\n",
    "\n",
    "-   `Hit`, to represent a single database hit. `Hit` objects are\n",
    "    contained within `QueryResult` and in each `QueryResult` there is\n",
    "    zero or more `Hit` objects.\n",
    "\n",
    "-   `HSP` (short for high-scoring pair), to represent region(s) of\n",
    "    significant alignments between query and hit sequences. `HSP`\n",
    "    objects are contained within `Hit` objects and each `Hit` has one or\n",
    "    more `HSP` objects.\n",
    "\n",
    "-   `HSPFragment`, to represent a single contiguous alignment between\n",
    "    query and hit sequences. `HSPFragment` objects are contained within\n",
    "    `HSP` objects. Most sequence search tools like BLAST and HMMER unify\n",
    "    `HSP` and `HSPFragment` objects as each `HSP` will only have a\n",
    "    single `HSPFragment`. However there are tools like BLAT and\n",
    "    Exonerate that produce `HSP` containing multiple `HSPFragment`.\n",
    "    Don’t worry if this seems a tad confusing now, we’ll elaborate more\n",
    "    on these two objects later on.\n",
    "\n",
    "These four objects are the ones you will interact with when you use\n",
    "`Bio.SearchIO`. They are created using one of the main `Bio.SearchIO`\n",
    "methods: `read`, `parse`, `index`, or `index_db`. The details of these\n",
    "methods are provided in later sections. For this section, we’ll only be\n",
    "using read and parse. These functions behave similarly to their\n",
    "`Bio.SeqIO` and `Bio.AlignIO` counterparts:\n",
    "\n",
    "-   `read` is used for search output files with a single query and\n",
    "    returns a `QueryResult` object\n",
    "\n",
    "-   `parse` is used for search output files with multiple queries and\n",
    "    returns a generator that yields `QueryResult` objects\n",
    "\n",
    "With that settled, let’s start probing each `Bio.SearchIO` object,\n",
    "beginning with `QueryResult`.\n",
    "\n",
    "### QueryResult {#sec:searchio-qresult}\n",
    "\n",
    "The QueryResult object represents a single search query and contains\n",
    "zero or more Hit objects. Let’s see what it looks like using the BLAST\n",
    "file we have:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blastn (2.3.0+)\n",
      "  Query: gi|8332116|gb|BE037100.1|BE037100 (1111)\n",
      "         MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold ac...\n",
      " Target: nt\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0      1  gi|731339628|ref|XM_010682658.1|  PREDICTED: Beta vulga...\n",
      "            1      1  gi|568824607|ref|XM_006466626.1|  PREDICTED: Citrus sin...\n",
      "            2      1  gi|568824605|ref|XM_006466625.1|  PREDICTED: Citrus sin...\n",
      "            3      1  gi|568824603|ref|XM_006466624.1|  PREDICTED: Citrus sin...\n",
      "            4      1  gi|568824601|ref|XM_006466623.1|  PREDICTED: Citrus sin...\n",
      "            5      1  gi|568824599|ref|XM_006466622.1|  PREDICTED: Citrus sin...\n",
      "            6      1  gi|567866318|ref|XM_006425719.1|  Citrus clementina hyp...\n",
      "            7      1  gi|567866316|ref|XM_006425718.1|  Citrus clementina hyp...\n",
      "            8      1  gi|567866314|ref|XM_006425717.1|  Citrus clementina hyp...\n",
      "            9      1  gi|567866312|ref|XM_006425716.1|  Citrus clementina hyp...\n",
      "           10      1  gi|590704208|ref|XM_007047033.1|  Theobroma cacao Cold-...\n",
      "           11      1  gi|590704205|ref|XM_007047032.1|  Theobroma cacao Cold-...\n",
      "           12      1  gi|694428700|ref|XM_009343631.1|  PREDICTED: Pyrus x br...\n",
      "           13      1  gi|694402986|ref|XM_009378191.1|  PREDICTED: Pyrus x br...\n",
      "           14      1  gi|743838297|ref|XM_011027373.1|  PREDICTED: Populus eu...\n",
      "           15      1  gi|743838293|ref|XM_011027372.1|  PREDICTED: Populus eu...\n",
      "           16      1  gi|595807351|ref|XM_007202530.1|  Prunus persica hypoth...\n",
      "           17      1  gi|566180892|ref|XM_006380679.1|  Populus trichocarpa c...\n",
      "           18      1  gi|764593175|ref|XM_004300526.2|  PREDICTED: Fragaria v...\n",
      "           19      1  gi|731440276|ref|XM_002274845.3|  PREDICTED: Vitis vini...\n",
      "           20      1  gi|349709091|emb|FQ378501.1|  Vitis vinifera clone SS0A...\n",
      "           21      1  gi|719997221|ref|XM_010256725.1|  PREDICTED: Nelumbo nu...\n",
      "           22      1  gi|645272858|ref|XM_008243375.1|  PREDICTED: Prunus mum...\n",
      "           23      1  gi|645272856|ref|XM_008243374.1|  PREDICTED: Prunus mum...\n",
      "           24      1  gi|848856318|ref|XM_013000712.1|  PREDICTED: Erythranth...\n",
      "           25      1  gi|255562758|ref|XM_002522339.1|  Ricinus communis COR4...\n",
      "           26      1  gi|658006068|ref|XM_008339966.1|  PREDICTED: Malus x do...\n",
      "           27      1  gi|802641059|ref|XM_012223735.1|  PREDICTED: Jatropha c...\n",
      "           28      1  gi|802641054|ref|XM_012223734.1|  PREDICTED: Jatropha c...\n",
      "           29      1  gi|823184330|ref|XM_012633708.1|  PREDICTED: Gossypium ...\n",
      "           ~~~\n",
      "           47      1  gi|729314350|ref|XM_010532905.1|  PREDICTED: Tarenaya h...\n",
      "           48      1  gi|731383573|ref|XM_002284686.2|  PREDICTED: Vitis vini...\n",
      "           49      1  gi|255762732|gb|GQ370517.1|  Salvia miltiorrhiza cold a...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/SearchIO/__init__.py:211: BiopythonExperimentalWarning: Bio.SearchIO is an experimental submodule which may undergo significant changes prior to its future official release.\n",
      "  BiopythonExperimentalWarning)\n"
     ]
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "blast_qresult = SearchIO.read('data/my_blast.xml', 'blast-xml')\n",
    "print(blast_qresult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We’ve just begun to scratch the surface of the object model, but you can\n",
    "see that there’s already some useful information. By invoking `print` on\n",
    "the `QueryResult` object, you can see:\n",
    "\n",
    "-   The program name and version (blastn version 2.2.27+)\n",
    "\n",
    "-   The query ID, description, and its sequence length (ID is 42291,\n",
    "    description is ‘mystery\\_seq’, and it is 61 nucleotides long)\n",
    "\n",
    "-   The target database to search against (refseq\\_rna)\n",
    "\n",
    "-   A quick overview of the resulting hits. For our query sequence,\n",
    "    there are 100 potential hits (numbered 0–99 in the table). For each\n",
    "    hit, we can also see how many HSPs it contains, its ID, and a\n",
    "    snippet of its description. Notice here that `Bio.SearchIO`\n",
    "    truncates the hit table overview, by showing only hits numbered\n",
    "    0–29, and then 97–99.\n",
    "\n",
    "Now let’s check our BLAT results using the same procedure as above:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blat (<unknown version>)\n",
      "  Query: mystery_seq (61)\n",
      "         <unknown description>\n",
      " Target: <unknown target>\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0     17  chr19  <unknown description>\n"
     ]
    }
   ],
   "source": [
    "blat_qresult = SearchIO.read('data/my_blat.psl', 'blat-psl')\n",
    "print(blat_qresult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You’ll immediately notice that there are some differences. Some of these\n",
    "are caused by the way PSL format stores its details, as you’ll see. The\n",
    "rest are caused by the genuine program and target database differences\n",
    "between our BLAST and BLAT searches:\n",
    "\n",
    "-   The program name and version. `Bio.SearchIO` knows that the program\n",
    "    is BLAT, but in the output file there is no information regarding\n",
    "    the program version so it defaults to ‘&lt;unknown version&gt;’.\n",
    "\n",
    "-   The query ID, description, and its sequence length. Notice here that\n",
    "    these details are slightly different from the ones we saw in BLAST.\n",
    "    The ID is ‘mystery\\_seq’ instead of 42991, there is no known\n",
    "    description, but the query length is still 61. This is actually a\n",
    "    difference introduced by the file formats themselves. BLAST\n",
    "    sometimes creates its own query IDs and uses your original ID as the\n",
    "    sequence description.\n",
    "\n",
    "-   The target database is not known, as it is not stated in the BLAT\n",
    "    output file.\n",
    "\n",
    "-   And finally, the list of hits we have is completely different. Here,\n",
    "    we see that our query sequence only hits the ‘chr19’ database entry,\n",
    "    but in it we see 17 HSP regions. This should not be surprising\n",
    "    however, given that we are using a different program, each with its\n",
    "    own target database.\n",
    "\n",
    "All the details you saw when invoking the `print` method can be accessed\n",
    "individually using Python’s object attribute access notation (a.k.a. the\n",
    "dot notation). There are also other format-specific attributes that you\n",
    "can access using the same method.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blastn 2.3.0+\n"
     ]
    }
   ],
   "source": [
    "print(\"%s %s\" % (blast_qresult.program, blast_qresult.version))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blat <unknown version>\n"
     ]
    }
   ],
   "source": [
    "print(\"%s %s\" % (blat_qresult.program, blat_qresult.version))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult.param_evalue_threshold    # blast-xml specific"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "For a complete list of accessible attributes, you can check each\n",
    "format-specific documentation. Here are the ones [for\n",
    "BLAST](http://biopython.org/DIST/docs/api/Bio.SearchIO.BlastIO-module.html)\n",
    "and for\n",
    "[BLAT](http://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html).\n",
    "\n",
    "Having looked at using `print` on `QueryResult` objects, let’s drill\n",
    "down deeper. What exactly is a `QueryResult`? In terms of Python\n",
    "objects, `QueryResult` is a hybrid between a list and a dictionary. In\n",
    "other words, it is a container object with all the convenient features\n",
    "of lists and dictionaries.\n",
    "\n",
    "Like Python lists and dictionaries, `QueryResult` objects are iterable.\n",
    "Each iteration returns a `Hit` object:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for hit in blast_qresult:\n",
    "    hit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To check how many items (hits) a `QueryResult` has, you can simply\n",
    "invoke Python’s `len` method:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(blast_qresult)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(blat_qresult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Like Python lists, you can retrieve items (hits) from a `QueryResult`\n",
    "using the slice notation:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Hit(id='gi|731339628|ref|XM_010682658.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult[0]        # retrieves the top hit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Hit(id='gi|255762732|gb|GQ370517.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult[-1]       # retrieves the last hit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To retrieve multiple hits, you can slice `QueryResult` objects using the\n",
    "slice notation as well. In this case, the slice will return a new\n",
    "`QueryResult` object containing only the sliced hits:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blastn (2.3.0+)\n",
      "  Query: gi|8332116|gb|BE037100.1|BE037100 (1111)\n",
      "         MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold ac...\n",
      " Target: nt\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0      1  gi|731339628|ref|XM_010682658.1|  PREDICTED: Beta vulga...\n",
      "            1      1  gi|568824607|ref|XM_006466626.1|  PREDICTED: Citrus sin...\n",
      "            2      1  gi|568824605|ref|XM_006466625.1|  PREDICTED: Citrus sin...\n"
     ]
    }
   ],
   "source": [
    "blast_slice = blast_qresult[:3]     # slices the first three hits\n",
    "print(blast_slice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Like Python dictionaries, you can also retrieve hits using the hit’s ID.\n",
    "This is particularly useful if you know a given hit ID exists within a\n",
    "search query results:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Hit(id='gi|731339628|ref|XM_010682658.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult['gi|731339628|ref|XM_010682658.1|']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You can also get a full list of `Hit` objects using `hits` and a full\n",
    "list of `Hit` IDs using `hit_keys`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Hit(id='gi|731339628|ref|XM_010682658.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|568824607|ref|XM_006466626.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|568824605|ref|XM_006466625.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|568824603|ref|XM_006466624.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|568824601|ref|XM_006466623.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|568824599|ref|XM_006466622.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|567866318|ref|XM_006425719.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|567866316|ref|XM_006425718.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|567866314|ref|XM_006425717.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|567866312|ref|XM_006425716.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|590704208|ref|XM_007047033.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|590704205|ref|XM_007047032.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|694428700|ref|XM_009343631.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|694402986|ref|XM_009378191.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|743838297|ref|XM_011027373.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|743838293|ref|XM_011027372.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|595807351|ref|XM_007202530.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|566180892|ref|XM_006380679.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|764593175|ref|XM_004300526.2|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|731440276|ref|XM_002274845.3|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|349709091|emb|FQ378501.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|719997221|ref|XM_010256725.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|645272858|ref|XM_008243375.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|645272856|ref|XM_008243374.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|848856318|ref|XM_013000712.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|255562758|ref|XM_002522339.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|658006068|ref|XM_008339966.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|802641059|ref|XM_012223735.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|802641054|ref|XM_012223734.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|823184330|ref|XM_012633708.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|697116717|ref|XM_009613983.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|698554552|ref|XM_009771948.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|565355903|ref|XM_006344753.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|659079297|ref|XM_008441961.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|659079295|ref|XM_008441960.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|828326432|ref|XM_004509143.2|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|729300104|ref|XM_010559563.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|922336892|ref|XM_013590982.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|955375440|ref|XM_003548859.3|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|778718146|ref|XM_011659511.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|778718144|ref|XM_004141977.2|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|955306249|ref|XM_003519866.3|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|747102191|ref|XM_011100949.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|723677217|ref|XM_004233368.2|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|703126872|ref|XM_010105390.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|658037046|ref|XM_008355863.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|225311746|dbj|AK326681.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|729314350|ref|XM_010532905.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|731383573|ref|XM_002284686.2|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps),\n",
       " Hit(id='gi|255762732|gb|GQ370517.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult.hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gi|731339628|ref|XM_010682658.1|',\n",
       " 'gi|568824607|ref|XM_006466626.1|',\n",
       " 'gi|568824605|ref|XM_006466625.1|',\n",
       " 'gi|568824603|ref|XM_006466624.1|',\n",
       " 'gi|568824601|ref|XM_006466623.1|',\n",
       " 'gi|568824599|ref|XM_006466622.1|',\n",
       " 'gi|567866318|ref|XM_006425719.1|',\n",
       " 'gi|567866316|ref|XM_006425718.1|',\n",
       " 'gi|567866314|ref|XM_006425717.1|',\n",
       " 'gi|567866312|ref|XM_006425716.1|',\n",
       " 'gi|590704208|ref|XM_007047033.1|',\n",
       " 'gi|590704205|ref|XM_007047032.1|',\n",
       " 'gi|694428700|ref|XM_009343631.1|',\n",
       " 'gi|694402986|ref|XM_009378191.1|',\n",
       " 'gi|743838297|ref|XM_011027373.1|',\n",
       " 'gi|743838293|ref|XM_011027372.1|',\n",
       " 'gi|595807351|ref|XM_007202530.1|',\n",
       " 'gi|566180892|ref|XM_006380679.1|',\n",
       " 'gi|764593175|ref|XM_004300526.2|',\n",
       " 'gi|731440276|ref|XM_002274845.3|',\n",
       " 'gi|349709091|emb|FQ378501.1|',\n",
       " 'gi|719997221|ref|XM_010256725.1|',\n",
       " 'gi|645272858|ref|XM_008243375.1|',\n",
       " 'gi|645272856|ref|XM_008243374.1|',\n",
       " 'gi|848856318|ref|XM_013000712.1|',\n",
       " 'gi|255562758|ref|XM_002522339.1|',\n",
       " 'gi|658006068|ref|XM_008339966.1|',\n",
       " 'gi|802641059|ref|XM_012223735.1|',\n",
       " 'gi|802641054|ref|XM_012223734.1|',\n",
       " 'gi|823184330|ref|XM_012633708.1|',\n",
       " 'gi|697116717|ref|XM_009613983.1|',\n",
       " 'gi|698554552|ref|XM_009771948.1|',\n",
       " 'gi|565355903|ref|XM_006344753.1|',\n",
       " 'gi|659079297|ref|XM_008441961.1|',\n",
       " 'gi|659079295|ref|XM_008441960.1|',\n",
       " 'gi|828326432|ref|XM_004509143.2|',\n",
       " 'gi|729300104|ref|XM_010559563.1|',\n",
       " 'gi|922336892|ref|XM_013590982.1|',\n",
       " 'gi|955375440|ref|XM_003548859.3|',\n",
       " 'gi|778718146|ref|XM_011659511.1|',\n",
       " 'gi|778718144|ref|XM_004141977.2|',\n",
       " 'gi|955306249|ref|XM_003519866.3|',\n",
       " 'gi|747102191|ref|XM_011100949.1|',\n",
       " 'gi|723677217|ref|XM_004233368.2|',\n",
       " 'gi|703126872|ref|XM_010105390.1|',\n",
       " 'gi|658037046|ref|XM_008355863.1|',\n",
       " 'gi|225311746|dbj|AK326681.1|',\n",
       " 'gi|729314350|ref|XM_010532905.1|',\n",
       " 'gi|731383573|ref|XM_002284686.2|',\n",
       " 'gi|255762732|gb|GQ370517.1|']"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult.hit_keys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "What if you just want to check whether a particular hit is present in\n",
    "the query results? You can do a simple Python membership test using the\n",
    "`in` keyword:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'gi|262205317|ref|NR_030195.1|' in blast_qresult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'gi|731339628|ref|XM_010682658.1|' in blast_qresult"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Sometimes, knowing whether a hit is present is not enough; you also want\n",
    "to know the rank of the hit. Here, the `index` method comes to the\n",
    "rescue:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_qresult.index('gi|595807351|ref|XM_007202530.1|')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Remember that we’re using Python’s indexing style here, which is\n",
    "zero-based. This means our hit above is ranked at no. 23, not 22.\n",
    "\n",
    "Also, note that the hit rank you see here is based on the native hit\n",
    "ordering present in the original search output file. Different search\n",
    "tools may order these hits based on different criteria.\n",
    "\n",
    "If the native hit ordering doesn’t suit your taste, you can use the\n",
    "`sort` method of the `QueryResult` object. It is very similar to\n",
    "Python’s `list.sort` method, with the addition of an option to create a\n",
    "new sorted `QueryResult` object or not.\n",
    "\n",
    "Here is an example of using `QueryResult.sort` to sort the hits based on\n",
    "each hit’s full sequence length. For this particular sort, we’ll set the\n",
    "`in_place` flag to `False` so that sorting will return a new\n",
    "`QueryResult` object and leave our initial object unsorted. We’ll also\n",
    "set the `reverse` flag to `True` so that we sort in descending order.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gi|731339628|ref|XM_010682658.1| 847\n",
      "gi|568824607|ref|XM_006466626.1| 878\n",
      "gi|568824605|ref|XM_006466625.1| 911\n",
      "gi|568824603|ref|XM_006466624.1| 894\n",
      "gi|568824601|ref|XM_006466623.1| 922\n"
     ]
    }
   ],
   "source": [
    "for hit in blast_qresult[:5]:   # id and sequence length of the first five hits\n",
    "    print(\"%s %i\" % (hit.id, hit.seq_len))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gi|955306249|ref|XM_003519866.3| 1388\n",
      "gi|659079295|ref|XM_008441960.1| 1245\n",
      "gi|567866316|ref|XM_006425718.1| 1202\n",
      "gi|566180892|ref|XM_006380679.1| 1182\n",
      "gi|565355903|ref|XM_006344753.1| 1153\n"
     ]
    }
   ],
   "source": [
    "sort_key = lambda hit: hit.seq_len\n",
    "sorted_qresult = blast_qresult.sort(key=sort_key, reverse=True, in_place=False)\n",
    "for hit in sorted_qresult[:5]:\n",
    "    print(\"%s %i\" % (hit.id, hit.seq_len))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The advantage of having the `in_place` flag here is that we’re\n",
    "preserving the native ordering, so we may use it again later. You should\n",
    "note that this is not the default behavior of `QueryResult.sort`,\n",
    "however, which is why we needed to set the `in_place` flag to `True`\n",
    "explicitly.\n",
    "\n",
    "At this point, you’ve known enough about `QueryResult` objects to make\n",
    "it work for you. But before we go on to the next object in the\n",
    "`Bio.SearchIO` model, let’s take a look at two more sets of methods that\n",
    "could make it even easier to work with `QueryResult` objects: the\n",
    "`filter` and `map` methods.\n",
    "\n",
    "If you’re familiar with Python’s list comprehensions, generator\n",
    "expressions or the built in `filter` and `map` functions, you’ll know\n",
    "how useful they are for working with list-like objects (if you’re not,\n",
    "check them out!). You can use these built in methods to manipulate\n",
    "`QueryResult` objects, but you’ll end up with regular Python lists and\n",
    "lose the ability to do more interesting manipulations.\n",
    "\n",
    "That’s why, `QueryResult` objects provide its own flavor of `filter` and\n",
    "`map` methods. Analogous to `filter`, there are `hit_filter` and\n",
    "`hsp_filter` methods. As their name implies, these methods filter its\n",
    "`QueryResult` object either on its `Hit` objects or `HSP` objects.\n",
    "Similarly, analogous to `map`, `QueryResult` objects also provide the\n",
    "`hit_map` and `hsp_map` methods. These methods apply a given function to\n",
    "all hits or HSPs in a `QueryResult` object, respectively.\n",
    "\n",
    "Let’s see these methods in action, beginning with `hit_filter`. This\n",
    "method accepts a callback function that checks whether a given `Hit`\n",
    "object passes the condition you set or not. In other words, the function\n",
    "must accept as its argument a single `Hit` object and returns `True` or\n",
    "`False`.\n",
    "\n",
    "Here is an example of using `hit_filter` to filter out `Hit` objects\n",
    "that only have one HSP:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filter_func = lambda hit: len(hit.hsps) > 1     # the callback function\n",
    "len(blast_qresult)      # no. of hits before filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filtered_qresult = blast_qresult.hit_filter(filter_func)\n",
    "len(filtered_qresult)   # no. of hits after filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for hit in filtered_qresult[:5]:    # quick check for the hit lengths\n",
    "    print(\"%s %i\" % (hit.id, len(hit.hsps)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "`hsp_filter` works the same as `hit_filter`, only instead of looking at\n",
    "the `Hit` objects, it performs filtering on the `HSP` objects in each\n",
    "hits.\n",
    "\n",
    "As for the `map` methods, they too accept a callback function as their\n",
    "arguments. However, instead of returning `True` or `False`, the callback\n",
    "function must return the modified `Hit` or `HSP` object (depending on\n",
    "whether you’re using `hit_map` or `hsp_map`).\n",
    "\n",
    "Let’s see an example where we’re using `hit_map` to rename the hit IDs:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XM_010682658.1\n",
      "XM_006466626.1\n",
      "XM_006466625.1\n",
      "XM_006466624.1\n",
      "XM_006466623.1\n"
     ]
    }
   ],
   "source": [
    "def map_func(hit):\n",
    "    hit.id = hit.id.split('|')[3]   # renames 'gi|301171322|ref|NR_035857.1|' to 'NR_035857.1'\n",
    "    return hit\n",
    "\n",
    "mapped_qresult = blast_qresult.hit_map(map_func)\n",
    "for hit in mapped_qresult[:5]:\n",
    "    print(hit.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Again, `hsp_map` works the same as `hit_map`, but on `HSP` objects\n",
    "instead of `Hit` objects.\n",
    "\n",
    "### Hit {#sec:searchio-hit}\n",
    "\n",
    "`Hit` objects represent all query results from a single database entry.\n",
    "They are the second-level container in the `Bio.SearchIO` object\n",
    "hierarchy. You’ve seen that they are contained by `QueryResult` objects,\n",
    "but they themselves contain `HSP` objects.\n",
    "\n",
    "Let’s see what they look like, beginning with our BLAST search:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Query: gi|8332116|gb|BE037100.1|BE037100\n",
      "       MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold accl...\n",
      "  Hit: gi|568824603|ref|XM_006466624.1| (894)\n",
      "       PREDICTED: Citrus sinensis cold-regulated 413 plasma membrane protein...\n",
      " HSPs: ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          #   E-value  Bit score    Span      Query range              Hit range\n",
      "       ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          0   6.9e-99     372.78     596         [63:655]              [108:697]\n"
     ]
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "blast_qresult = SearchIO.read('data/my_blast.xml', 'blast-xml')\n",
    "blast_hit = blast_qresult[3]    # fourth hit from the query result\n",
    "print(blast_hit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You see that we’ve got the essentials covered here:\n",
    "\n",
    "-   The query ID and description is present. A hit is always tied to a\n",
    "    query, so we want to keep track of the originating query as well.\n",
    "    These values can be accessed from a hit using the `query_id` and\n",
    "    `query_description` attributes.\n",
    "\n",
    "-   We also have the unique hit ID, description, and full\n",
    "    sequence lengths. They can be accessed using `id`, `description`,\n",
    "    and `seq_len`, respectively.\n",
    "\n",
    "-   Finally, there’s a table containing quick information about the HSPs\n",
    "    this hit contains. In each row, we’ve got the important HSP details\n",
    "    listed: the HSP index, its e-value, its bit score, its span (the\n",
    "    alignment length including gaps), its query coordinates, and its\n",
    "    hit coordinates.\n",
    "\n",
    "Now let’s contrast this with the BLAT search. Remember that in the BLAT\n",
    "search we had one hit with 17 HSPs.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Query: mystery_seq\n",
      "       <unknown description>\n",
      "  Hit: chr19 (59128983)\n",
      "       <unknown description>\n",
      " HSPs: ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          #   E-value  Bit score    Span      Query range              Hit range\n",
      "       ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          0         ?          ?       ?           [0:61]    [54204480:54204541]\n",
      "          1         ?          ?       ?           [0:61]    [54233104:54264463]\n",
      "          2         ?          ?       ?           [0:61]    [54254477:54260071]\n",
      "          3         ?          ?       ?           [1:61]    [54210720:54210780]\n",
      "          4         ?          ?       ?           [0:60]    [54198476:54198536]\n",
      "          5         ?          ?       ?           [0:61]    [54265610:54265671]\n",
      "          6         ?          ?       ?           [0:61]    [54238143:54240175]\n",
      "          7         ?          ?       ?           [0:60]    [54189735:54189795]\n",
      "          8         ?          ?       ?           [0:61]    [54185425:54185486]\n",
      "          9         ?          ?       ?           [0:60]    [54197657:54197717]\n",
      "         10         ?          ?       ?           [0:61]    [54255662:54255723]\n",
      "         11         ?          ?       ?           [0:61]    [54201651:54201712]\n",
      "         12         ?          ?       ?           [8:60]    [54206009:54206061]\n",
      "         13         ?          ?       ?          [10:61]    [54178987:54179038]\n",
      "         14         ?          ?       ?           [8:61]    [54212018:54212071]\n",
      "         15         ?          ?       ?           [8:51]    [54234278:54234321]\n",
      "         16         ?          ?       ?           [8:61]    [54238143:54238196]\n"
     ]
    }
   ],
   "source": [
    "blat_qresult = SearchIO.read('data/my_blat.psl', 'blat-psl')\n",
    "blat_hit = blat_qresult[0]      # the only hit\n",
    "print(blat_hit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Here, we’ve got a similar level of detail as with the BLAST hit we saw\n",
    "earlier. There are some differences worth explaining, though:\n",
    "\n",
    "-   The e-value and bit score column values. As BLAT HSPs do not have\n",
    "    e-values and bit scores, the display defaults to ‘?’.\n",
    "\n",
    "-   What about the span column? The span values is meant to display the\n",
    "    complete alignment length, which consists of all residues and any\n",
    "    gaps that may be present. The PSL format do not have this\n",
    "    information readily available and `Bio.SearchIO` does not attempt to\n",
    "    try guess what it is, so we get a ‘?’ similar to the e-value and bit\n",
    "    score columns.\n",
    "\n",
    "In terms of Python objects, `Hit` behaves almost the same as Python\n",
    "lists, but contain `HSP` objects exclusively. If you’re familiar with\n",
    "lists, you should encounter no difficulties working with the `Hit`\n",
    "object.\n",
    "\n",
    "Just like Python lists, `Hit` objects are iterable, and each iteration\n",
    "returns one `HSP` object it contains:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for hsp in blast_hit:\n",
    "    hsp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You can invoke `len` on a `Hit` to see how many `HSP` objects it has:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(blast_hit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(blat_hit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You can use the slice notation on `Hit` objects, whether to retrieve\n",
    "single `HSP` or multiple `HSP` objects. Like `QueryResult`, if you slice\n",
    "for multiple `HSP`, a new `Hit` object will be returned containing only\n",
    "the sliced `HSP` objects:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "HSP(hit_id='chr19', query_id='mystery_seq', 1 fragments)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hit[0]                 # retrieve single items"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sliced_hit = blat_hit[4:9]  # retrieve multiple items\n",
    "len(sliced_hit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Query: mystery_seq\n",
      "       <unknown description>\n",
      "  Hit: chr19 (59128983)\n",
      "       <unknown description>\n",
      " HSPs: ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          #   E-value  Bit score    Span      Query range              Hit range\n",
      "       ----  --------  ---------  ------  ---------------  ---------------------\n",
      "          0         ?          ?       ?           [0:60]    [54198476:54198536]\n",
      "          1         ?          ?       ?           [0:61]    [54265610:54265671]\n",
      "          2         ?          ?       ?           [0:61]    [54238143:54240175]\n",
      "          3         ?          ?       ?           [0:60]    [54189735:54189795]\n",
      "          4         ?          ?       ?           [0:61]    [54185425:54185486]\n"
     ]
    }
   ],
   "source": [
    "print(sliced_hit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You can also sort the `HSP` inside a `Hit`, using the exact same\n",
    "arguments like the sort method you saw in the `QueryResult` object.\n",
    "\n",
    "Finally, there are also the `filter` and `map` methods you can use on\n",
    "`Hit` objects. Unlike in the `QueryResult` object, `Hit` objects only\n",
    "have one variant of `filter` (`Hit.filter`) and one variant of `map`\n",
    "(`Hit.map`). Both of `Hit.filter` and `Hit.map` work on the `HSP`\n",
    "objects a `Hit` has.\n",
    "\n",
    "### HSP {#sec:searchio-hsp}\n",
    "\n",
    "`HSP` (high-scoring pair) represents region(s) in the hit sequence that\n",
    "contains significant alignment(s) to the query sequence. It contains the\n",
    "actual match between your query sequence and a database entry. As this\n",
    "match is determined by the sequence search tool’s algorithms, the `HSP`\n",
    "object contains the bulk of the statistics computed by the search tool.\n",
    "This also makes the distinction between `HSP` objects from different\n",
    "search tools more apparent compared to the differences you’ve seen in\n",
    "`QueryResult` or `Hit` objects.\n",
    "\n",
    "Let’s see some examples from our BLAST and BLAT searches. We’ll look at\n",
    "the BLAST HSP first:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: gi|8332116|gb|BE037100.1|BE037100 MP14H09 MP Mesembryanthemum cr...\n",
      "        Hit: gi|731339628|ref|XM_010682658.1| PREDICTED: Beta vulgaris subsp....\n",
      "Query range: [86:679] (1)\n",
      "  Hit range: [80:677] (1)\n",
      "Quick stats: evalue 1.2e-108; bitscore 405.24\n",
      "  Fragments: 1 (597 columns)\n",
      "     Query - TTGGCCATGAAAACTGATCAATTGGCCGTGGCTAATATGATCGATTCCGATATCAATGA~~~TGTAG\n",
      "             ||||||||||||||||| ||| ||||  |||||||| |||| ||||  ||||| |||||~~~|||||\n",
      "       Hit - TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGATTATGATATGAATGA~~~TGTAG\n"
     ]
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "blast_qresult = SearchIO.read('data/my_blast.xml', 'blast-xml')\n",
    "blast_hsp = blast_qresult[0][0]    # first hit, first hsp\n",
    "print(blast_hsp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Just like `QueryResult` and `Hit`, invoking `print` on an `HSP` shows\n",
    "its general details:\n",
    "\n",
    "-   There are the query and hit IDs and descriptions. We need these to\n",
    "    identify our `HSP`.\n",
    "\n",
    "-   We’ve also got the matching range of the query and hit sequences.\n",
    "    The slice notation we’re using here is an indication that the range\n",
    "    is displayed using Python’s indexing style (zero-based, half open).\n",
    "    The number inside the parenthesis denotes the strand. In this case,\n",
    "    both sequences have the plus strand.\n",
    "\n",
    "-   Some quick statistics are available: the e-value and bitscore.\n",
    "\n",
    "-   There is information about the HSP fragments. Ignore this for now;\n",
    "    it will be explained later on.\n",
    "\n",
    "-   And finally, we have the query and hit sequence alignment itself.\n",
    "\n",
    "These details can be accessed on their own using the dot notation, just\n",
    "like in `QueryResult` and `Hit`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(86, 679)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.query_range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1684e-108"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.evalue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "They’re not the only attributes available, though. `HSP` objects come\n",
    "with a default set of properties that makes it easy to probe their\n",
    "various details. Here are some examples:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.hit_start         # start coordinate of the hit sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def map_func(hit):\n",
    "    hit.id = hit.id.split('|')[3]   # renames 'gi|301171322|ref|NR_035857.1|' to 'NR_035857.1'\n",
    "    return hitblast_hsp.query_span        # how many residues in the query sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "597"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.aln_span          # how long the alignment is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Check out the `HSP`\n",
    "[documentation](http://biopython.org/DIST/docs/api/Bio.SearchIO._model.hsp-module.html)\n",
    "for a full list of these predefined properties.\n",
    "\n",
    "Furthermore, each sequence search tool usually computes its own\n",
    "statistics / details for its `HSP` objects. For example, an XML BLAST\n",
    "search also outputs the number of gaps and identical residues. These\n",
    "attributes can be accessed like so:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.gap_num       # number of gaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "449"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.ident_num     # number of identical residues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "These details are format-specific; they may not be present in other\n",
    "formats. To see which details are available for a given sequence search\n",
    "tool, you should check the format’s documentation in `Bio.SearchIO`.\n",
    "Alternatively, you may also use `.__dict__.keys()` for a quick list of\n",
    "what’s available:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['_items', 'pos_num', 'gap_num', 'bitscore', 'bitscore_raw', 'ident_num', 'evalue'])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.__dict__.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Finally, you may have noticed that the `query` and `hit` attributes of\n",
    "our HSP are not just regular strings:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('TTGGCCATGAAAACTGATCAATTGGCCGTGGCTAATATGATCGATTCCGATATC...TAG', DNAAlphabet()), id='gi|8332116|gb|BE037100.1|BE037100', name='aligned query sequence', description=\"MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold acclimation protein, mRNA sequence\", dbxrefs=[])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGATTATGATATG...TAG', DNAAlphabet()), id='gi|731339628|ref|XM_010682658.1|', name='aligned hit sequence', description='PREDICTED: Beta vulgaris subsp. vulgaris cold-regulated 413 plasma membrane protein 2 (LOC104895996), mRNA', dbxrefs=[])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_hsp.hit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "They are `SeqRecord` objects you saw earlier in\n",
    "Section \\[chapter:SeqRecord\\]! This means that you can do all sorts of\n",
    "interesting things you can do with `SeqRecord` objects on `HSP.query`\n",
    "and/or `HSP.hit`.\n",
    "\n",
    "It should not surprise you now that the `HSP` object has an `alignment`\n",
    "property which is a `MultipleSeqAlignment` object:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DNAAlphabet() alignment with 2 rows and 597 columns\n",
      "TTGGCCATGAAAACTGATCAATTGGCCGTGGCTAATATGATCGA...TAG gi|8332116|gb|BE037100.1|BE037100\n",
      "TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGA...TAG gi|731339628|ref|XM_010682658.1|\n"
     ]
    }
   ],
   "source": [
    "print(blast_hsp.aln)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Having probed the BLAST HSP, let’s now take a look at HSPs from our BLAT\n",
    "results for a different kind of HSP. As usual, we’ll begin by invoking\n",
    "`print` on it:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: mystery_seq <unknown description>\n",
      "        Hit: chr19 <unknown description>\n",
      "Query range: [0:61] (1)\n",
      "  Hit range: [54204480:54204541] (1)\n",
      "Quick stats: evalue ?; bitscore ?\n",
      "  Fragments: 1 (? columns)\n"
     ]
    }
   ],
   "source": [
    "blat_qresult = SearchIO.read('data/my_blat.psl', 'blat-psl')\n",
    "blat_hsp = blat_qresult[0][0]       # first hit, first hsp\n",
    "print(blat_hsp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Some of the outputs you may have already guessed. We have the query and\n",
    "hit IDs and descriptions and the sequence coordinates. Values for evalue\n",
    "and bitscore is ‘?’ as BLAT HSPs do not have these attributes. But The\n",
    "biggest difference here is that you don’t see any sequence alignments\n",
    "displayed. If you look closer, PSL formats themselves do not have any\n",
    "hit or query sequences, so `Bio.SearchIO` won’t create any sequence or\n",
    "alignment objects. What happens if you try to access `HSP.query`,\n",
    "`HSP.hit`, or `HSP.aln`? You’ll get the default values for these\n",
    "attributes, which is `None`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.hit is None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.query is None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.aln is None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This does not affect other attributes, though. For example, you can\n",
    "still access the length of the query or hit alignment. Despite not\n",
    "displaying any attributes, the PSL format still have this information so\n",
    "`Bio.SearchIO` can extract them:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.query_span     # length of query match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.hit_span       # length of hit match"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Other format-specific attributes are still present as well:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.score          # PSL score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.mismatch_num   # the mismatch column"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "So far so good? Things get more interesting when you look at another\n",
    "‘variant’ of HSP present in our BLAT results. You might recall that in\n",
    "BLAT searches, sometimes we get our results separated into ‘blocks’.\n",
    "These blocks are essentially alignment fragments that may have some\n",
    "intervening sequence between them.\n",
    "\n",
    "Let’s take a look at a BLAT HSP that contains multiple blocks to see how\n",
    "`Bio.SearchIO` deals with this:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: mystery_seq <unknown description>\n",
      "        Hit: chr19 <unknown description>\n",
      "Query range: [0:61] (1)\n",
      "  Hit range: [54233104:54264463] (1)\n",
      "Quick stats: evalue ?; bitscore ?\n",
      "  Fragments: ---  --------------  ----------------------  ----------------------\n",
      "               #            Span             Query range               Hit range\n",
      "             ---  --------------  ----------------------  ----------------------\n",
      "               0               ?                  [0:18]     [54233104:54233122]\n",
      "               1               ?                 [18:61]     [54264420:54264463]\n"
     ]
    }
   ],
   "source": [
    "blat_hsp2 = blat_qresult[0][1]      # first hit, second hsp\n",
    "print(blat_hsp2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "What’s happening here? We still some essential details covered: the IDs\n",
    "and descriptions, the coordinates, and the quick statistics are similar\n",
    "to what you’ve seen before. But the fragments detail is all different.\n",
    "Instead of showing ‘Fragments: 1’, we now have a table with two data\n",
    "rows.\n",
    "\n",
    "This is how `Bio.SearchIO` deals with HSPs having multiple fragments. As\n",
    "mentioned before, an HSP alignment may be separated by intervening\n",
    "sequences into fragments. The intervening sequences are not part of the\n",
    "query-hit match, so they should not be considered part of query nor hit\n",
    "sequence. However, they do affect how we deal with sequence coordinates,\n",
    "so we can’t ignore them.\n",
    "\n",
    "Take a look at the hit coordinate of the HSP above. In the `Hit range:`\n",
    "field, we see that the coordinate is `[54233104:54264463]`. But looking\n",
    "at the table rows, we see that not the entire region spanned by this\n",
    "coordinate matches our query. Specifically, the intervening region spans\n",
    "from `54233122` to `54264420`.\n",
    "\n",
    "Why then, is the query coordinates seem to be contiguous, you ask? This\n",
    "is perfectly fine. In this case it means that the query match is\n",
    "contiguous (no intervening regions), while the hit match is not.\n",
    "\n",
    "All these attributes are accessible from the HSP directly, by the way:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(54233104, 54264463)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_range         # hit start and end coordinates of the entire HSP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(54233104, 54233122), (54264420, 54264463)]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_range_all     # hit start and end coordinates of each fragment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31359"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_span          # hit span of the entire HSP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[18, 43]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_span_all      # hit span of each fragment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(54233122, 54264420)]"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_inter_ranges  # start and end coordinates of intervening regions in the hit sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[31298]"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.hit_inter_spans   # span of intervening regions in the hit sequence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Most of these attributes are not readily available from the PSL file we\n",
    "have, but `Bio.SearchIO` calculates them for you on the fly when you\n",
    "parse the PSL file. All it needs are the start and end coordinates of\n",
    "each fragment.\n",
    "\n",
    "What about the `query`, `hit`, and `aln` attributes? If the HSP has\n",
    "multiple fragments, you won’t be able to use these attributes as they\n",
    "only fetch single `SeqRecord` or `MultipleSeqAlignment` objects.\n",
    "However, you can use their `*_all` counterparts: `query_all`, `hit_all`,\n",
    "and `aln_all`. These properties will return a list containing\n",
    "`SeqRecord` or `MultipleSeqAlignment` objects from each of the HSP\n",
    "fragment. There are other attributes that behave similarly, i.e. they\n",
    "only work for HSPs with one fragment. Check out the `HSP`\n",
    "[documentation](http://biopython.org/DIST/docs/api/Bio.SearchIO._model.hsp-module.html)\n",
    "for a full list.\n",
    "\n",
    "Finally, to check whether you have multiple fragments or not, you can\n",
    "use the `is_fragmented` property like so:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp2.is_fragmented     # BLAT HSP with 2 fragments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blat_hsp.is_fragmented      # BLAT HSP from earlier, with one fragment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Before we move on, you should also know that we can use the slice\n",
    "notation on `HSP` objects, just like `QueryResult` or `Hit` objects.\n",
    "When you use this notation, you’ll get an `HSPFragment` object in\n",
    "return, the last component of the object model.\n",
    "\n",
    "### HSPFragment {#sec:searchio-hspfragment}\n",
    "\n",
    "`HSPFragment` represents a single, contiguous match between the query\n",
    "and hit sequences. You could consider it the core of the object model\n",
    "and search result, since it is the presence of these fragments that\n",
    "determine whether your search have results or not.\n",
    "\n",
    "In most cases, you don’t have to deal with `HSPFragment` objects\n",
    "directly since not that many sequence search tools fragment their HSPs.\n",
    "When you do have to deal with them, what you should remember is that\n",
    "`HSPFragment` objects were written with to be as compact as possible. In\n",
    "most cases, they only contain attributes directly related to sequences:\n",
    "strands, reading frames, alphabets, coordinates, the sequences\n",
    "themselves, and their IDs and descriptions.\n",
    "\n",
    "These attributes are readily shown when you invoke `print` on an\n",
    "`HSPFragment`. Here’s an example, taken from our BLAST search:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: gi|8332116|gb|BE037100.1|BE037100 MP14H09 MP Mesembryanthemum cr...\n",
      "        Hit: gi|731339628|ref|XM_010682658.1| PREDICTED: Beta vulgaris subsp....\n",
      "Query range: [86:679] (1)\n",
      "  Hit range: [80:677] (1)\n",
      "  Fragments: 1 (597 columns)\n",
      "     Query - TTGGCCATGAAAACTGATCAATTGGCCGTGGCTAATATGATCGATTCCGATATCAATGA~~~TGTAG\n",
      "             ||||||||||||||||| ||| ||||  |||||||| |||| ||||  ||||| |||||~~~|||||\n",
      "       Hit - TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGATTATGATATGAATGA~~~TGTAG\n"
     ]
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "blast_qresult = SearchIO.read('data/my_blast.xml', 'blast-xml')\n",
    "blast_frag = blast_qresult[0][0][0]    # first hit, first hsp, first fragment\n",
    "print(blast_frag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "At this level, the BLAT fragment looks quite similar to the BLAST\n",
    "fragment, save for the query and hit sequences which are not present:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: mystery_seq <unknown description>\n",
      "        Hit: chr19 <unknown description>\n",
      "Query range: [0:61] (1)\n",
      "  Hit range: [54204480:54204541] (1)\n",
      "  Fragments: 1 (? columns)\n"
     ]
    }
   ],
   "source": [
    "blat_qresult = SearchIO.read('data/my_blat.psl', 'blat-psl')\n",
    "blat_frag = blat_qresult[0][0][0]    # first hit, first hsp, first fragment\n",
    "print(blat_frag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In all cases, these attributes are accessible using our favorite dot\n",
    "notation. Some examples:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "86"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_frag.query_start      # query start coordinate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_frag.hit_strand       # hit sequence strand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGATTATGATATG...TAG', DNAAlphabet()), id='gi|731339628|ref|XM_010682658.1|', name='aligned hit sequence', description='PREDICTED: Beta vulgaris subsp. vulgaris cold-regulated 413 plasma membrane protein 2 (LOC104895996), mRNA', dbxrefs=[])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blast_frag.hit              # hit sequence, as a SeqRecord object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "A note about standards and conventions {#sec:searchio-standards}\n",
    "--------------------------------------\n",
    "\n",
    "Before we move on to the main functions, there is something you ought to\n",
    "know about the standards `Bio.SearchIO` uses. If you’ve worked with\n",
    "multiple sequence search tools, you might have had to deal with the many\n",
    "different ways each program deals with things like sequence coordinates.\n",
    "It might not have been a pleasant experience as these search tools\n",
    "usually have their own standards. For example, one tools might use\n",
    "one-based coordinates, while the other uses zero-based coordinates. Or,\n",
    "one program might reverse the start and end coordinates if the strand is\n",
    "minus, while others don’t. In short, these often creates unnecessary\n",
    "mess must be dealt with.\n",
    "\n",
    "We realize this problem ourselves and we intend to address it in\n",
    "`Bio.SearchIO`. After all, one of the goals of `Bio.SearchIO` is to\n",
    "create a common, easy to use interface to deal with various search\n",
    "output files. This means creating standards that extend beyond the\n",
    "object model you just saw.\n",
    "\n",
    "Now, you might complain, “Not another standard!”. Well, eventually we\n",
    "have to choose one convention or the other, so this is necessary. Plus,\n",
    "we’re not creating something entirely new here; just adopting a standard\n",
    "we think is best for a Python programmer (it is Biopython, after all).\n",
    "\n",
    "There are three implicit standards that you can expect when working with\n",
    "`Bio.SearchIO`:\n",
    "\n",
    "-   The first one pertains to sequence coordinates. In `Bio.SearchIO`,\n",
    "    all sequence coordinates follows Python’s coordinate style:\n",
    "    zero-based and half open. For example, if in a BLAST XML output file\n",
    "    the start and end coordinates of an HSP are 10 and 28, they would\n",
    "    become 9 and 28 in `Bio.SearchIO`. The start coordinate becomes 9\n",
    "    because Python indices start from zero, while the end coordinate\n",
    "    remains 28 as Python slices omit the last item in an interval.\n",
    "\n",
    "-   The second is on sequence coordinate orders. In `Bio.SearchIO`,\n",
    "    start coordinates are always less than or equal to end coordinates.\n",
    "    This isn’t always the case with all sequence search tools, as some\n",
    "    of them have larger start coordinates when the sequence strand\n",
    "    is minus.\n",
    "\n",
    "-   The last one is on strand and reading frame values. For strands,\n",
    "    there are only four valid choices: `1` (plus strand), `-1` (minus\n",
    "    strand), `0` (protein sequences), and `None` (no strand). For\n",
    "    reading frames, the valid choices are integers from `-3` to `3` and\n",
    "    `None`.\n",
    "\n",
    "Note that these standards only exist in `Bio.SearchIO` objects. If you\n",
    "write `Bio.SearchIO` objects into an output format, `Bio.SearchIO` will\n",
    "use the format’s standard for the output. It does not force its standard\n",
    "over to your output file.\n",
    "\n",
    "Reading search output files {#sec:searchio-input}\n",
    "---------------------------\n",
    "\n",
    "There are two functions you can use for reading search output files into\n",
    "`Bio.SearchIO` objects: `read` and `parse`. They’re essentially similar\n",
    "to `read` and `parse` functions in other submodules like `Bio.SeqIO` or\n",
    "`Bio.AlignIO`. In both cases, you need to supply the search output file\n",
    "name and the file format name, both as Python strings. You can check the\n",
    "documentation for a list of format names `Bio.SearchIO` recognizes.\n",
    "\n",
    "`Bio.SearchIO.read` is used for reading search output files with only\n",
    "one query and returns a `QueryResult` object. You’ve seen `read` used in\n",
    "our previous examples. What you haven’t seen is that `read` may also\n",
    "accept additional keyword arguments, depending on the file format.\n",
    "\n",
    "Here are some examples. In the first one, we use `read` just like\n",
    "previously to read a BLAST tabular output file. In the second one, we\n",
    "use a keyword argument to modify so it parses the BLAST tabular variant\n",
    "with comments in it:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "qresult = SearchIO.read('data/tab_2226_tblastn_003.txt', 'blast-tab')\n",
    "qresult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qresult2 = SearchIO.read('data/tab_2226_tblastn_007.txt', 'blast-tab', comments=True)\n",
    "qresult2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "These keyword arguments differs among file formats. Check the format\n",
    "documentation to see if it has keyword arguments that modifies its\n",
    "parser’s behavior.\n",
    "\n",
    "As for the `Bio.SearchIO.parse`, it is used for reading search output\n",
    "files with any number of queries. The function returns a generator\n",
    "object that yields a `QueryResult` object in each iteration. Like\n",
    "`Bio.SearchIO.read`, it also accepts format-specific keyword arguments:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gi|16080617|ref|NP_391444.1|\n",
      "gi|11464971:4-101\n"
     ]
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "qresults = SearchIO.parse('data/tab_2226_tblastn_001.txt', 'blast-tab')\n",
    "for qresult in qresults:\n",
    "    print(qresult.id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "random_s00\n",
      "gi|16080617|ref|NP_391444.1|\n",
      "gi|11464971:4-101\n"
     ]
    }
   ],
   "source": [
    "qresults2 = SearchIO.parse('data/tab_2226_tblastn_005.txt', 'blast-tab', comments=True)\n",
    "for qresult in qresults2:\n",
    "    print(qresult.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Dealing with large search output files with indexing {#sec:searchio-index}\n",
    "----------------------------------------------------\n",
    "\n",
    "Sometimes, you’re handed a search output file containing hundreds or\n",
    "thousands of queries that you need to parse. You can of course use\n",
    "`Bio.SearchIO.parse` for this file, but that would be grossly\n",
    "inefficient if you need to access only a few of the queries. This is\n",
    "because `parse` will parse all queries it sees before it fetches your\n",
    "query of interest.\n",
    "\n",
    "In this case, the ideal choice would be to index the file using\n",
    "`Bio.SearchIO.index` or `Bio.SearchIO.index_db`. If the names sound\n",
    "familiar, it’s because you’ve seen them before in\n",
    "Section \\[sec:SeqIO-index\\]. These functions also behave similarly to\n",
    "their `Bio.SeqIO` counterparts, with the addition of format-specific\n",
    "keyword arguments.\n",
    "\n",
    "Here are some examples. You can use `index` with just the filename and\n",
    "format name:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|']"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "idx = SearchIO.index('data/tab_2226_tblastn_001.txt', 'blast-tab')\n",
    "sorted(idx.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx['gi|16080617|ref|NP_391444.1|']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idx.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Or also with the format-specific keyword argument:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|', 'random_s00']"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx = SearchIO.index('data/tab_2226_tblastn_005.txt', 'blast-tab', comments=True)\n",
    "sorted(idx.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx['gi|16080617|ref|NP_391444.1|']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idx.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Or with the `key_function` argument, as in `Bio.SeqIO`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GI|11464971:4-101', 'GI|16080617|REF|NP_391444.1|']"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "key_function = lambda id: id.upper()    # capitalizes the keys\n",
    "idx = SearchIO.index('data/tab_2226_tblastn_001.txt', 'blast-tab', key_function=key_function)\n",
    "sorted(idx.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx['GI|16080617|REF|NP_391444.1|']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idx.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "`Bio.SearchIO.index_db` works like as `index`, only it writes the query\n",
    "offsets into an SQLite database file.\n",
    "\n",
    "Writing and converting search output files {#sec:searchio-write}\n",
    "------------------------------------------\n",
    "\n",
    "It is occasionally useful to be able to manipulate search results from\n",
    "an output file and write it again to a new file. `Bio.SearchIO` provides\n",
    "a `write` function that lets you do exactly this. It takes as its\n",
    "arguments an iterable returning `QueryResult` objects, the output\n",
    "filename to write to, the format name to write to, and optionally some\n",
    "format-specific keyword arguments. It returns a four-item tuple, which\n",
    "denotes the number or `QueryResult`, `Hit`, `HSP`, and `HSPFragment`\n",
    "objects that were written.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 239, 277, 277)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "qresults = SearchIO.parse('data/mirna.xml', 'blast-xml')     # read XML file\n",
    "SearchIO.write(qresults, 'results.tab', 'blast-tab')    # write to tabular file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You should note different file formats require different attributes of\n",
    "the `QueryResult`, `Hit`, `HSP` and `HSPFragment` objects. If these\n",
    "attributes are not present, writing won’t work. In other words, you\n",
    "can’t always write to the output format that you want. For example, if\n",
    "you read a BLAST XML file, you wouldn’t be able to write the results to\n",
    "a PSL file as PSL files require attributes not calculated by BLAST (e.g.\n",
    "the number of repeat matches). You can always set these attributes\n",
    "manually, if you really want to write to PSL, though.\n",
    "\n",
    "Like `read`, `parse`, `index`, and `index_db`, `write` also accepts\n",
    "format-specific keyword arguments. Check out the documentation for a\n",
    "complete list of formats `Bio.SearchIO` can write to and their\n",
    "arguments.\n",
    "\n",
    "Finally, `Bio.SearchIO` also provides a `convert` function, which is\n",
    "simply a shortcut for `Bio.SearchIO.parse` and `Bio.SearchIO.write`.\n",
    "Using the convert function, our example above would be:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 239, 277, 277)"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SearchIO\n",
    "SearchIO.convert('data/mirna.xml', 'blast-xml', 'results.tab', 'blast-tab')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
